Analyse spatiale et territoriale
Formation Carthageo-Geoprisme 2023
Syllabus
Objectifs de la semaine
- Développer les connaissances en analyse quantitative des données dans l’espace
- Exploiter des données individuelles (demandes de valeurs foncières, recensement de la population)
- Approfondir l’utilisation du logiciel R
- Travailler en groupe et préparer une restitution commune
- Répondre à une “commande” dans un temps limité
Organisation de la semaine
Planning des intervenants et contenu des journées
- 09/10 (matin, Pierre Pistre) : Organisation de la semaine & présentation des données
- 09/10 (après-midi, Claude Grasland) : Prise en main des données de recensement (fichier LOGEMT)
- 10/10 (matin, Pierre Pistre) : Prise en main des données de valeurs foncières (DVF)
- 10/10 (après-midi, Claude Grasland) : extractions des données pour la commande, et approfondissements statistiques & cartographiques
- 11/10 : Accompagnement des projets de groupe (Hadrien Commenges)
- 12/10 : Accompagnement des projets de groupe (Hadrien Commenges)
- 13/10 : Finalisation des projets de groupe (matin) et restitution orale (après-midi)
Horaires et planning des salles (Campus des Grands Moulins, Université Paris Cité)
- Horaires indicatifs : 9h30-12h30 et 13h30-17h.
- Salles : du 09/10 au 11/10 : salle de cours 209 (Bâtiment Olympe de Gouges, 2ième étage) ; 12/10 : salle informatique 375 (Bâtiment Olympe de Gouges, 3ième étage) ; 13/10 : salle de cours 209 (Bâtiment Olympe de Gouges, 2ième étage)
Données utilisées
(sources) Base de données principales :
- Fichier “Logements oridinaires” (LOGEMT) 2019 (localisations : IRIS, communes…), produit par l’INSEE à partir du recensement de la population. Données en accés libre : https://www.insee.fr/fr/statistiques/6544344?sommaire=6456104
- “Demandes de valeurs foncières” (version Etalab et extraction avec Opendatasoft), produit par le Ministère de l’Économie, des Finances et de la Souveraineté industrielle et numérique : https://www.data.gouv.fr/fr/datasets/demandes-de-valeurs-foncieres-geolocalisees/
(sources) Fichiers complémentaires :
- Table d’appartenance géographique des communes (INSEE) : https://www.insee.fr/fr/information/2028028
- Shapefile ADMIN-EXPRESS-COG® (IGN) des découpage administratifs de la France métropolitaine : https://geoservices.ign.fr/adminexpress#telechargement
- Shapefile Contours IRIS® (IGN) : https://geoservices.ign.fr/contoursiris
Exercice d’application en groupe
Périmètres de l’exercice
- Objet d’étude : les dynamiques immobilières dans les principales agglomérations administratives de la “mégarégion parisienne” (hors Paris)
- Population d’étude : les logements vendus récemment et leurs habitants
- Espace(s) d’étude : les principales agglomérations ayant notamment le statut de “métrople” ou de “communauté urbaine” (hors Paris)
- Echelles et mailles géographiques : département, commune, IRIS
Consignes de la “commande”
- Contexte général : il est connu que Paris influence la dynamique des territoires bien au-delà de ses territoires les plus proches, à commencer par les principales agglomérations urbaines de la “mégarégion parisienne” (au sens élargi de : https://atlas-paris-mega-region.univ-rouen.fr/). Le secteur immobilier est particulièrement concerné du fait d’une progression sensible des prix des logements dans la plupart de ces territoires au cours à minima de la dernière décennie et d’un potentiel renforcement avec la crise COVID par une progression des installations résidentielles. Pour autant, la géographie des dynamiques immobilières internes à ces agglomérations urbaines reste mal connue de même que les facteurs explicatifs (influence parisienne, spécificités locales et régionales…) et les tendances récentes.
- Commande : après avoir fait le constat de dynamiques immobilières souvent proches, les principales agglomérations urbaines de la “mégarégion parisienne” ont décidé de s’associer pour demander aux différentes antennes régionales de l’Insee la réalisation d’études rigoureuses et fines spatialement sur la situation immobilière de leurs différents territoires de compétence et périphéries proches. Le bilan de ces études sera présenté dans des publications synthétiques, intelligibles par le plus grand nombre et par les acteurs concernés, sous la forme de plusieurs 4 pages - comme produits régulièrement par l’Insee pour publier les résultats de ses études (par exemple, file:///home/pierre/T%C3%A9l%C3%A9chargements/ip1715.pdf ; file:///home/pierre/T%C3%A9l%C3%A9chargements/IR134_Notaires-IPLA_1T-23.pdf).
Modalités de l’exercice, restitution du travail et organisation durant la semaine
- Constitution de 6 (voire 7) groupes de 4-5 étudiant.e.s, mélangeant les profils de Master (= Lundi matin)
- Réflexion sur un angle problématique pour chaque cas d’étude, ainsi que les possibilités d’analyse et d’exploration des données : variables pertinentes, traitements envisagés… (= Lundi à Mercredi)
- Réalisation des analyses statistiques et cartographiques (= Mardi à Jeudi)
- Organisation, mise en forme des analyses et rédaction de la note écrite (= Mardi à Vendredi matin)
- Présentation orale (environ 10 minutes, sans autre support visuel que la note écrite) devant un “jury” composé des intervenants de la semaine (= Vendredi après-midi)
Cas d’étude
- Caen (communauté urbaine)
- Le Mans (communauté urbaine)
- Orléans (métropole)
- Reims (communauté urbaine)
- Rouen (métropole)
- Tours (métropole)
En option supplémentaire : Amiens (communauté d’agglomération)
Format du rendu (type 4 pages Insee)
- Introduction et définition (thématique et statistique) de l’objet : thème et cas d’étude
- Etat des lieux général des prix immobiliers dans le cas d’étude (données : DVF ; échelle : “agglomération” (rayon 40km) ; maille : commune)
- Approfondissements sur les prix immobiliers par type de biens ou de localisation (données : DVF ; échelles : “agglomération” (rayon 40km) et zooms ; mailles : commune, IRIS)
- Profils des acheteurs au sens des nouveaux arrivants (moins de 3 ans) qui sont propriétaires (données : recensement de la population ; échelles : “agglomération” (rayon 40km) et zooms ; mailles : commune, IRIS)
RP1 : Accquisition
Introduction
Ce programme permet de préparer l’ensemble des données statistiques et cartographiques utiles pour l’étude d’une agglomération urbaine à partir du fichier ménage du recensement de population 2019. Il charge de très gros fichiers contenus dans un dossier temporaire tmp et en extrait les informations utiles pour les placer dans un dossier data/nomagg . Il a besoin en entrée du code et du nom de la commune centre ainsi que du chemin vers le dossier de stockage des résultats
Etape 1 : Données administratives
On commence par récupérer les informations adminstratives sur la zone d’étude
Chargement du fichier administatif
Identification des codes de l’aire urbaine et de l’EPCI
On extrait la ligne correspondant à notre commune centre et on en déduit le code de son epci et de son aire urbaine
Extraction des informations
On extrait toutes les donnes relatives à notre aire urbaine et on les sauvegarde.
Etape 2 : Données géométriques
On va maintenant extraire les données géométriques relatives au contour des communes et des IRIS
Communes
On charge le fichier des communes de France et on extrait uniquement celles de la zone d’étude.
## Communes de France
mapcom <- st_read("tmp/maps/COMMUNE.shp", quiet=T)
## Aire urbaine
AA_adm<-readRDS(paste0(myrep,"AA_adm.RDS"))
AA_mapcom <- mapcom %>% filter(INSEE_COM %in% AA_adm$CODGEO)
saveRDS(AA_mapcom,paste0(myrep,"AA_mapcom.RDS"))
## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_mapcom <- mapcom %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_mapcom,paste0(myrep,"EPCI_mapcom.RDS"))
## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_mapcom <- mapcom %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_mapcom,paste0(myrep,"CTR_mapcom.RDS"))AA_mapcom<-readRDS(paste0(myrep,"AA_mapcom.RDS"))
EPCI_mapcom<-readRDS(paste0(myrep,"EPCI_mapcom.RDS"))
CTR_mapcom<-readRDS(paste0(myrep,"CTR_mapcom.RDS"))
plot(AA_mapcom$geometry, col="lightyellow", main="zone d'étude")
plot(EPCI_mapcom$geometry,col="orange", add=T)
plot(CTR_mapcom$geometry,col="red", add=T)IRIS
On ne procède à l’extraction des IRIS que pour l’EPCI et la commune centre
## Communes de France
mapiris <- st_read("tmp/maps/CONTOURS-IRIS.shp", quiet=T)
## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_mapiris <- mapiris %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_mapiris,paste0(myrep,"EPCI_mapiris.RDS"))
## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_mapiris <- mapiris %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_mapiris,paste0(myrep,"CTR_mapiris.RDS"))EPCI_mapiris<-readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
CTR_mapiris<-readRDS(paste0(myrep,"CTR_mapiris.RDS"))
plot(EPCI_mapiris$geometry, col="lightyellow", main="EPCI",border="gray50", lwd=0.4)
plot(CTR_mapiris$geometry,col="orange", add=T, border="gray50", lwd=0.4)
plot(EPCI_mapcom$geometry,col=NA, border="black",lwd=1,add=T)Etape 3 : logements ordinaires en 2019
Nous partirons des fichiers détail de l’INSEE car, à la différence des tableaux prédéfinis, ils permettent virtuellement toutes les formes de croisement d’indicateurs. Ils sont évidemment beaucoup plus volumineux, mais ce sera justement l’occasion pour les étudiants en data mining d’être confrontés à des problèmes d’optimisation et de big data. On trouve leur description détaillée sur le site de l’INSEE
Compte-tenu de la taille des fichiers, nous travaillerons sur une partition de la France en zones
- Zone A : Région Île-de-France (région 11) ;
- Zone B : Régions Centre-Val de Loire (région 24), Bourgogne-Franche-Comté (région 27), Normandie (région 28) et Hauts-de-France (région 32) ;
- Zone C : Régions Grand Est (région 44), Pays de la Loire (région 52) et Bretagne (région 53);
- Zone D : Régions Nouvelle-Aquitaine (région 75) et Occitanie (région 76) ;
- Zone E : Régions Auvergne-Rhônes-Alpes (région 84), Provence-Alpes-Côte d’Azur (région 93), Corse (région 94), Guadeloupe (région 01), Martinique (région 02), Guyane (région 03) et La Réunion (région 04).
Lecture du fichier des ménages
Nous commençons par lire le fichiers de la zone qui nous intéresse (au format .csv). Nous utilisons pour cela la fonction fread du package data.table qui est très rapide mais qui commet une erreur sur le code communal qu’on doit corriger
Extraction des ménages de la zone d’étude
## Aire urbaine
AA_adm<-readRDS(paste0(myrep,"AA_adm.RDS"))
AA_men <- men %>% filter(INSEE_COM %in% AA_adm$CODGEO)
saveRDS(AA_men,paste0(myrep,"AA_men.RDS"))
## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_men <- men %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_men,paste0(myrep,"EPCI_men.RDS"))
## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_men <- men %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_men,paste0(myrep,"CTR_men.RDS"))Chargement des méta-données
Il ne faut surtout pas oublier le fichier des métadonnées qui va permettre de recoder facilement tous les facteurs et de décoder les chiffres correspondant aux classes. On va donc le transformer au format R puis l’enregistrer également dans le dossier data. On élimine toutefois les codes des variables COMMUNES, IRIS et TRIRIS qui prennent trop de place.
RP2 : Agrégation
Le format sf (spatial features)
La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois
- un tableau de données (l’équivalent du fichier .dbf)
- une géométrie (l’équivalent du fichier .shp)
- une projection (l’équivalent du fichier .prj)
Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou
dans d’autres formats standards comme GeoJson, la première tâche
consiste donc à les convertir au formt sf afin de pouvoir les utiliser
facilement dans R. L’importation se fait à l’aide de l’instruction
st_read en indiquant juste le nom du fichier .shp à
charger. Les autres fichiers (.dbf ou .proj) seront lus également et
intégrés dans l’objet qui hérite de la double classe data.frame
et sf.
Etapes de préparation des données
Dans notre exemple, nous allons suivre les étapes suivantes :
- Préparer les données statistiques par IRIS dans un data.frame
- Charger un fonds de carte par IRIS au format sf
- Effectuer une jointure entre les deux fichiers par le code IRIS
- Sauvegarder le résultat
- Agréger les données statistiques et géométriques par commune
- Sauvegarder le résultat.
Préparer les données statistiques
On importe le fichier des individus et on corrige le code IRIS des communes non divisées en IRIS.
myrep<-"data/Rennes/"
tab_ind<-readRDS(paste0(myrep,"EPCI_men.RDS"))
tab_ind$IRIS[tab_ind$IRIS=="ZZZZZZZZZ"]<-
paste0(tab_ind$INSEE_COM[tab_ind$IRIS=="ZZZZZZZZZ"],"0000")
head(tab_ind[,1:5],3)FALSE COMMUNE ARM IRIS ACHL AEMM
FALSE 1: 35001 ZZZZZ 350010101 B12 1976
FALSE 2: 35001 ZZZZZ 350010101 B12 1977
FALSE 3: 35001 ZZZZZ 350010101 B12 0
Selectionner les lignes et colonnes
On décide de ne sélectionner que les ménages de propriétaires installés depuis moins de 3 ans et on retient comme variable le type de logement que l’on recode en deux catégories
tab_ind2 <- tab_ind %>% filter(as.numeric(ANEM) < 3, # durée d'occupation < 3 ans
STOCD =="10", # propriétaire de son logement
TYPL %in% c(1,2)) %>% # Maison ou Appartement
mutate(TYPL = as.factor(TYPL)) %>%
select(IPONDL, IRIS, TYPL)
levels(tab_ind2$TYPL)<-c("Maison","Appt")
knitr::kable(head(tab_ind2,5),digits=2)| IPONDL | IRIS | TYPL |
|---|---|---|
| 1.07 | 350010101 | Maison |
| 1.03 | 350010101 | Maison |
| 1.07 | 350010101 | Maison |
| 0.99 | 350010101 | Maison |
| 1.11 | 350010101 | Maison |
Agréger les données
On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :
tab_long<- tab_ind2 %>%
group_by(IRIS,TYPL)%>%
summarise(NB=sum(IPONDL))
knitr::kable(head(tab_long,5),digits=2)| IRIS | TYPL | NB |
|---|---|---|
| 350010101 | Maison | 36.13 |
| 350010101 | Appt | 9.71 |
| 350010102 | Maison | 106.53 |
| 350010102 | Appt | 57.60 |
| 350220000 | Maison | 21.89 |
Pivoter le tableau
Puis on fait “pivoter” le tableau pour l’obtenir en format large :
tab_large <- tab_long %>% pivot_wider(id_cols = IRIS,
names_from = TYPL,
names_prefix = "T_",
values_from = NB,
values_fill = 0)
knitr::kable(head(tab_large,5),digits=2)| IRIS | T_Maison | T_Appt |
|---|---|---|
| 350010101 | 36.13 | 9.71 |
| 350010102 | 106.53 | 57.60 |
| 350220000 | 21.89 | 0.99 |
| 350240101 | 142.85 | 41.34 |
| 350240102 | 34.58 | 10.50 |
Ajouter de nouvelles variables
On ajoute de nouvelles variables telles que le nombre total de nouveaux propriétaires et la part des propriétaires de maisons parmi eux.
tab<- tab_large %>% mutate(T_TOT = T_Maison+T_Appt,
Maison_pct = 100*T_Maison/T_TOT)
knitr::kable(head(tab,5),digits=c(0,0,0,0,1))| IRIS | T_Maison | T_Appt | T_TOT | Maison_pct |
|---|---|---|---|---|
| 350010101 | 36 | 10 | 46 | 78.8 |
| 350010102 | 107 | 58 | 164 | 64.9 |
| 350220000 | 22 | 1 | 23 | 95.7 |
| 350240101 | 143 | 41 | 184 | 77.6 |
| 350240102 | 35 | 10 | 45 | 76.7 |
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de nouveux propriétaires en maison individuelle.
programme
p <- ggplot(tab) + aes (x = Maison_pct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90, 100)) +
scale_x_continuous("% de Maisons") +
scale_y_continuous("Nombre d'IRIS") +
ggtitle(label = "Type de logement des nouveaux propriétaires",
subtitle = "Source : INSEE, RP 2019 - Métropole de Rennes")
pCharger les données géométriques
On importe le fichier des iris de l’agglomération qui est de type sf (spatial features)
map_iris <- readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
map_iris<-map_iris[,c(5,6,3,2)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")
class(map_iris)FALSE [1] "sf" "data.frame"
| IRIS | NOM_IRIS | COM | NOM_COM | geometry |
|---|---|---|---|---|
| 352381204 | Les Cloteaux | Rennes | 35238 | MULTIPOLYGON (((350568 6786… |
| 352750000 | Saint-Gilles | Saint-Gilles | 35275 | MULTIPOLYGON (((337255.9 67… |
Visualisation du fonds iris avec sf
On peut facilement faire un plot de la colonne geometry du fichier sf
Jointure des données IRIS et du fonds de carte
programme
map_iris_tab<-merge(map_iris,tab,
by.x="IRIS",by.y="IRIS",
all.x=T,all.y=F)
knitr::kable(head(map_iris_tab,3),digits=2)| IRIS | NOM_IRIS | COM | NOM_COM | T_Maison | T_Appt | T_TOT | Maison_pct | geometry |
|---|---|---|---|---|---|---|---|---|
| 350010101 | Ouest | Acigné | 35001 | 36.13 | 9.71 | 45.84 | 78.82 | MULTIPOLYGON (((360989.4 67… |
| 350010102 | Est | Acigné | 35001 | 106.53 | 57.60 | 164.13 | 64.90 | MULTIPOLYGON (((364373.3 67… |
| 350220000 | Bécherel | Bécherel | 35022 | 21.89 | 0.99 | 22.88 | 95.65 | MULTIPOLYGON (((333530.4 68… |
Cartographie rapide avec sf
Une astuce pour visualiser rapidement une carte choropléthe avec le seul packahe sf :
Sauvegarde du fichier par IRIS
Agrégation statistique + géométriques
Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”
Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.
Agrégation des IRIS en communes
L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries
Agrégation des iris en communes
résultat statistique
| COM | NOM_COM | T_Maison | T_Appt | T_TOT | Maison_pct |
|---|---|---|---|---|---|
| Acigné | 35001 | 143 | 67 | 210 | 67.9 |
| Betton | 35024 | 291 | 82 | 373 | 78.0 |
| Bourgbarré | 35032 | 162 | 24 | 186 | 87.2 |
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de nouveaux propriétaires en maison par commune.
p <- ggplot(map_com_tab) +aes (x = Maison_pct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90, 100)) +
scale_x_continuous("% de Maisons") +
scale_y_continuous("Nombre de communes") +
ggtitle(label = "Type de logement des nouveaux propriétaires",
subtitle = "Source : INSEE, RP 2019, Agglomération de Rennes")
p RP3 : Cartographie
Le package map_sf
Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.
On trouvera la documentation du package mapsf à l’adresse suivante :
Création d’un template cartographique
Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :
tracé d’un fonds de carte vierge
La fonction mf_map() avec le paramètre
type = "base"permet de tracer une carte vide
Superposition de couches
On peut toutefois ajouter toute une série de paramètres
supplémentaire (col=, border=,
lwd=, …) et superposer plusieurs fonds de carte avec le
paramètre add = TRUE. L’ajout de la fonction
layout permet de rajouter un cadre une légende.
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray50",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=0.6,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude",
credits = "Sources : IGN et INSEE")Ajout d’un thème
On peut finalement modifier l’ensemble de la carte en lui ajoutant
une instruction mf_theme() qui peut reprendre des styles
existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”,
“candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”,
“barcelona”) mais aussi créer ses propres thèmes
#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray50",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=0.6,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude",
credits = "Sources : IGN et INSEE")Ajout de texte
On peut ajouter une couche de texte avec la fonction
mf_label(). Par exemple, on va ajouter à la carte
précédente le nom des communes
programme
mf_theme("agolalight")
# Trace les Iris avec des paramètres
mf_map(map_iris,
type = "base",
col = "lightyellow",
border="gray50",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="red",
lwd=0.6,
add = TRUE)
# Ajoute une couche de labels
mf_label(map_com,
var="COM",
cex=0.4,
col="blue",
overlap = FALSE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude",
credits = "Sources : IGN et INSEE")Carte de stock
Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.
Dans le package mapsf, on réalise ce type de carte à
l’aide de la fonction mf_map()en lui donnant le paramètre
type="prop".
On va tenter à titre d’exemple de représenter la distribution du nombre de nouveaux propriétaires.
Carte de stock minimale
Les instructions minimales sont les suivantes :
# Trace les contours des communes
mf_map(x= map_iris,
type = "base")
# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris,
type ="prop",
var = "T_TOT",
add=TRUE)Mais le résultat est peu satisfaisant car les cercles sont trop
grands. Il faut en pratique toujours effectuer un réglage de ceux-ci
avec l’instruction inches=
Carte de stock habillée
mf_theme("agolalight")
mf_map(map_iris, type = "base",
col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="black",lwd=1,add = TRUE)
mf_map(map_iris, var = "T_TOT",type = "prop",
inches = 0.05, col = "red",leg_pos = "left",
leg_title = "Nombre de ménages", add=TRUE)
mf_layout(title = "Distribution des nouveaux propriétaires",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte choroplèthe
Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.
La fonction du package mapsf adaptée aux variables
d’intensité est la fonction mf_map()munie du paramètre
type = "choro".
On va prendre l’exemple du % de maisons parmi les nouveaux propriétaires
Carte choroplèthe minimale
Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).
Carte choroplèthe habillée
On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.
mybreaks = c(0, 10,20,30,40,50,
60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5),
pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris, var = "Maison_pct",
type = "choro",
breaks = mybreaks,pal = mypal,
border=NA,
col_na = "gray80",leg_title = "% maisons",
leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com, type = "base", col = NA,
border="black",lwd=1,add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte stock + choroplèthe (1)
On peut combiner les deux modes cartographiques par superposition :
mf_theme("agolalight")
# Choisit les classes
mybreaks = c(0,20,40,60,80,100)
# Trace la carte choroplèthe
mf_map(
x = map_iris,
var = "Maison_pct",
breaks = mybreaks,
# pal=mypal,
type = "choro",
border="white",
col_na = "gray80",
lwd=0.3,
leg_title = "% Maisons",
leg_val_rnd = 0,
)
# Ajoute les cercles proportionnels
mf_map(
x =map_iris,
var = "T_Maison",
type = "prop",
inches = 0.06,
col = "red",
leg_pos = "right",
leg_title = "Effectif",
add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="black",
lwd=1,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les nouveaux propriétaires de maison",
frame = TRUE,
credits = "Sources : IGN et INSEE")Carte stock + choroplèthe (2)
Mais les cercles dissimuent alors les plages de couleur, aussi on
peut utiliser le type prop_choro qui place la variable
choroplèthe à l’intérieur des cercles
mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens","Reds"))
mf_map(map_iris, type = "base",
col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris, var = c("T_TOT", "Maison_pct"),
inches = 0.06, col_na = "grey", pal=mypal,
breaks = mybreaks, nbreaks = 4, lwd = 0.1,
leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
leg_title = c("Nouveaux propriétaires", "% Maison"),
add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons",
frame = TRUE, credits = "Sources : IGN et INSEE")DVF1 : Acquisition
Introduction
Objectif
On se propose de collecter l’ensemble des ventes de maison ou d’appartement dans un rayon de 50 km autour d’une ville. On choisit ici Rennes en prenant la place du Parlement de Bretagne comme centre de référence
Données dvf
Les données sur les dvf sont disponibles sur le site publicopendatasof:
https://public.opendatasoft.com
On s’intéresse plus particulièrement aux données dvf géolocalisées accessibles ici :
Sélection des informations
Le fichier comporte 28 millions d’enregistrement ce qui est évidemment beaucoup … On va donc procéder à une sélection en se servant des différents onglets disponibles.
A titre d’exemple, nous allons essayer de télécharger les enregistrements vérifiant les conditions suivantes :
- localisation dans un rayon de 60 km autour d’Amiens
- type de transaction : Vente
- type de bien : Maison ou appartement
Nous constatons qu’il y a 133294 enregistrements vérifiant ces conditions.
Lien de téléchargement
Nous pouvons récupérer les données soit au format .csv
pour les analyse statistiques, soit au format .geojson pour
les analyses spatiales. Mais dans ce cas on ne peut pas dépasser la
taille de 50 000 enregistrement. Il vaut donc mieux télécharger au
format .csv puisque ce fichier comporte les latitudes et longitudes ce
qui nous permettra de créer un fichier cartographique a posteriori
On peut effectuer le téléchargement depuis le site web ou bien juste enregistrer le lien de téléchargement et effectuer l’opération dans R.
Pour trouver le lien on passe la souris au dessus du lien “télécharger les données au format .csv / seulement les enregistrements sélectionnées” et on effectue un click droit suivi de “copier le lien”
Construction d’une API
On récupère le lien de téléchargement du fichier .csv :
A première vue ceci paraît assez obscur mais on réalise assez vite que cette URL reprend l’ensemble des spécifications envoyées pour sélectionner nos enregistrements. On peut souligner en gras les parties utiles :
- &refine=nature_mutation%3A%22Vente%22&
- &refine=type_local%3A%22Maison%22
- &refine=type_local%3A%22Appartement%22
- &where=(distance(%60geo_point%60%2C%20geom%27
- POINT(2.26968%2049.90807)%27%2C%2060037m%20))
Nous allons donc créer un lien qui extrait automatiquement les dvf en modifiant juste la position du point central et le rayon. On prend cette fois-ci nos paramètres de Rennes
myurl<-paste0("https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/exports/csv?lang=fr&refine=nature_mutation%3A%22Vente%22&refine=type_local%3A%22Maison%22&refine=type_local%3A%22Appartement%22&facet=facet(name%3D%22nature_mutation%22%2C%20disjunctive%3Dtrue)&facet=facet(name%3D%22type_local%22%2C%20disjunctive%3Dtrue)&where=(distance(%60geo_point%60%2C%20geom%27POINT(", latctr,"%20",lonctr,")%27%2C%20",rayon,"m))&timezone=Europe%2FParis&use_labels=true&delimiter=%3B")Téléchargement
On lance la requête avec l’instruction download.file()
et on stocke le résultat dans notre répertoire de travail :
essai de l’URL ‘https://public.opendatasoft.com/…’ downloaded 38.4 MB
L’opération s’execute en un peu moins d’une minute dans cet exemple. Mais il peut arriver qu’elle échoue lorsque le serveur est saturé et dans ce cas il faut recommencer.
DVF2 : Nettoyage
Introduction
Objectif
On se propose de nettoyer les donénes relative à l’ensemble des ventes de maison ou d’appartement dans un rayon de 50 km autour d’une ville. On rappelle les paramètres qui ont été utilisés pour acquérir ces données.
Paramètres
Merci à Boris et Florent !
Nous allons effectuer le nettoyage en nous appuyant sur les propositions de Boris Mericskay et Florent Demoraes parues dans https://journals.openedition.org/cybergeo/39583
Les auteurs fournissent en effet un accès à tous les programmes utilisés dans leur article via une archive github : https://github.com/ESO-Rennes/Analyse-Donnees-DVF
Chargement du jeu de données
On adapte le programme de Boris et Florent en utilisant la fonction
fread() du package data.table car elle est
plus rapide/
Etape 1 : Sélection des mutations de type “Ventes” de “Maison” et “Appartement’
Programme
Remarque
Etape inutile si on a déjà fait ce double filtrage à l’aide de l’API
Etape 2 : Sélection et renommage des variables
etape2 <- etape1bis %>%
select(id = `Identifiant de mutation (Etalab)`,
disposition = `Numéro de disposition` ,
parcelle= `Identifiant de la parcelle cadastrale`,
date = `Date de la mutation`,
nature = `Nature de la mutation`,
INSEE_COM = `Code INSEE de la commune`,
INSEE_DEP = `Code INSEE du département` ,
type = `Type de local` ,
surface = `Surface réelle du bâti` ,
piece = `Nombre de pièces principales`,
prix = `Valeur foncière` ,
latitude = `Latitude`,
longitude = `Longitude`)Etape 3 : Remplacement des cellules vides par des NA et suppression des NA
Etape 4 : Sélections des mutations simples
Regrouper les transactions selon l’ID, la surface et la vente
Etape 5 : Jointure attributaire pour récupérer les informations de la mutation
jointure
Modification des formats des colonnes
Suppression des valeurs aberrantes
Fixer un seuil minimal et maximal des prix (percentiles)
1% 99%
17000 616320
Fixer un seuil minimal et maximal des prix (histogramme)
Créer deux jeux de données (maisons / appartements)
Appartement Maison
42003 68734
Etape 6 : Sélection des bornes de prix et de surface
Etape 8 : Sélection des bornes de prix au m2
Fixer un seuil minimal des prix au m2 (percentile)
1% 99%
387.3057 5825.8048
Etape 9 : Corrections diverses
Transformation de la date en année
Arrondir les variables numériques
Etape 10 : Exportation
Mise en ordre des lignes et colonnes
Ecrire le jeu de données final en csv
DVF3 : Cartographie
Retour sur sf
Nous revenons sur le package sf (spatial features)
que nous avons déjà rencontré au moment de la création de cartes
thématiques par IRIS ou communes à l’aide du package
mapsf.
Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des ventes foncières au dessus du fonds de carte des IRIS ou communes.
Mais il pourra aussi servir de base à des cartographies
dynamiques permettant de placer les points correspondant aux
ventes sur des réseaux de rue et plus généralement sur des “tuiles”
cartographiques permettant d’effectur des zoom. On utilisera à cet effet
d’autres packages comme leaflet ou sa version simplifiée
mapview.
Zone d’étude
On définit les paramètres de la zone d’étude en ajoutant le lien vers les dvf nettoyées
Choix de communes
On choisit plus précisément de travailler sur les communes de l’agglmomération
map_com <- readRDS(paste0(myrep,"EPCI_mapcom.RDS")) %>%
select(INSEE_COM, NOM, geometry)
plot(map_com$geometry,
col="lightyellow")Sélection des dvf nettoyées
Nous reprenons le fichier nettoyé de localisation des ventes de maison et d’appartement établi au chapitre précédent et nous ne conservons que 9 variables ainsi que les ventes situées dans les communes de la zone choisie
dvf <- readRDS(paste0(myrep,mydvf)) %>%
select(annee, INSEE_COM,
type, prix, surface, prixm2,
latitude, longitude) %>%
filter(INSEE_COM %in% map_com$INSEE_COM)
kable(head(dvf,3))| annee | INSEE_COM | type | prix | surface | prixm2 | latitude | longitude |
|---|---|---|---|---|---|---|---|
| 2014 | 35238 | Appartement | 81000 | 66 | 1227 | 48.11526 | -1.692243 |
| 2014 | 35238 | Appartement | 98500 | 33 | 2985 | 48.11689 | -1.681832 |
| 2014 | 35238 | Appartement | 121000 | 62 | 1952 | 48.10974 | -1.705489 |
Transformation des dvf en fichier sf
On ajoute une colonne geometry en utilisant la fonction
st_as_sf()du package sf. On attribue à la
projection le code EPSG = 4326 qui correspond à un fichier latitude
longitude.
map_dvf <- st_as_sf(dvf, coords = c("longitude","latitude"))
st_crs(map_dvf)<- 4326
kable(head(map_dvf,3))| annee | INSEE_COM | type | prix | surface | prixm2 | geometry |
|---|---|---|---|---|---|---|
| 2014 | 35238 | Appartement | 81000 | 66 | 1227 | POINT (-1.692243 48.11526) |
| 2014 | 35238 | Appartement | 98500 | 33 | 2985 | POINT (-1.681832 48.11689) |
| 2014 | 35238 | Appartement | 121000 | 62 | 1952 | POINT (-1.705489 48.10974) |
Test de superposition
On tente de superposer les deux cartes des ventes et des communes.
Commune puis dvf
Vérification de la projection
Nous savons que les coordonnées du fichier DVF sont en latitude longitude (EPSG=4326). Mais la projection de notre fonds des communes est en Lambert-93 (EPSG = 2154)
Coordinate Reference System:
User input: RGF93 v1 / Lambert-93
wkt:
PROJCRS["RGF93 v1 / Lambert-93",
BASEGEOGCRS["RGF93 v1",
DATUM["Reseau Geodesique Francais 1993 v1",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4171]],
CONVERSION["Lambert-93",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",46.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",3,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",49,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",44,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",700000,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",6600000,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["France - onshore and offshore, mainland and Corsica."],
BBOX[41.15,-9.86,51.56,10.38]],
ID["EPSG",2154]]
A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS
Ajustement des deux projections
On transforme donc le fichier des DVF en Lambert-93 pour obtenir une adéquation avec le fichier des communes.
Et cette fois-ci on réussit à superposer les deux fonds de carte
Cartographie ponctuelle avec mapsf
On peut utiliser les fonctions de mapsf vues précédemment. Par exemple repérer les ventes de maison et d’appartement :
mf_map(map_com, type="base",
col="lightyellow")
mf_map(map_dvf, type = "typo",
var= "type",add =T,
border=NA,
cex=0.3,
leg_title = "Bien vendu")
mf_label(map_com,
var="NOM",cex = 0.4,
overlap = F,
col="black",
halo=T)
mf_layout(title = "Ventes de maisons et d'appartement 2013-2020",
credits = "Sources : IGN, INSEE, DVF",
scale = T)Cartographie ponctuelle avec mapview
Le package mapview pemet de créer facilement des cartes interactives ù l’on pourra visualiser les ventes sur un fonds de carte de son choix (OSRM, ESRI, …). On peut l’utiliser par exemple pour voirune commune précise
On peut améliorer sensiblement la visualisation à l’aide de paramètres présentés sur le site web du package https://r-spatial.github.io/mapview/
library(mapview)
# Choix de la zone d'étude
map_com_sel<-map_com %>% filter(INSEE_COM == "35196")
map_dvf_sel<-map_dvf %>% filter(INSEE_COM == "35196")
# Choix des tuiles
mapviewOptions(basemaps = c("OpenStreetMap" ,
"Esri.WorldImagery"))
# Première couche
map1 <- mapview(map_com_sel,
col.regions="lightyellow",
alpha.regions = 0.4)
# Deuxième couche
map2<-mapview(map_dvf_sel,
zcol = "type",
cex = 4,
col.regions=c("blue","red"),
alpha.regions = 0.4)
# Assemblage
map1+map2Agrégation par communes
On peut évidemment agréger les informations sur les valeurs foncières par communes. Par exemple compter le nombre de maisons vendues et leur prix moyen au m2
dvf_by_com <- dvf %>%
group_by(INSEE_COM) %>%
filter(type == "Maison") %>%
summarise(nb = n(),
med_prixm2 = median(prixm2))
kable(head(dvf_by_com,3))| INSEE_COM | nb | med_prixm2 |
|---|---|---|
| 35001 | 400 | 2417.0 |
| 35022 | 66 | 1369.5 |
| 35024 | 745 | 2576.0 |
On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :
On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :
mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_com_dvf, type = "choro",
var = "med_prixm2",
breaks=mybreaks, pal = mypal,
border="gray", lwd=0.3,
leg_title = "prix en €/m2", leg_val_rnd = 0,
leg_pos = "left")
mf_map( map_com_dvf, type = "prop",
var = "nb", inches = 0.12,
col="black",border="white",
leg_title = "nb. de ventes",
leg_pos = "topleft")
mf_layout(title = "Les ventes de maison de 2013 à 2020",
frame = TRUE, credits = "Sources : IGN et DVF",
arrow = F )Cartographie par IRIS
Supposons maintenant que l’on souhaite cartographier les prix des maisons par IRIS à l’intérieur de la commune de Rennes. Comment faire dans la mesure où l’IRIS n’est pas mentionné dans le fichier DVF ? Réponse : en superposant les deux fonds de carte et en les intersectant.
Superposition
map_iris_ctr<-readRDS(paste0(myrep,"CTR_mapiris.RDS")) %>%
select("NOM_COM","INSEE_COM","CODE_IRIS","NOM_IRIS","geometry")
map_dvf_ctr <- map_dvf %>% filter(INSEE_COM %in% map_iris_ctr$INSEE_COM)
plot(map_iris_ctr$geometry, col="lightyellow")
plot(map_dvf_ctr$geometry,col="red",pch=21,cex=0.1,add=T)Intersection
On peut faire l’opération avec un SIG ou bien avec le package
sf en utilisant sa fonction st_intersect()
qui va ajouter tous les attributs du fichier IRIS au fichier DVF :
| annee | INSEE_COM | type | prix | surface | prixm2 | NOM_COM | INSEE_COM.1 | CODE_IRIS | NOM_IRIS | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 62 | 2014 | 35238 | Appartement | 80000 | 44 | 1818 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350454.3 6786503) |
| 130 | 2014 | 35238 | Appartement | 70000 | 55 | 1273 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350613.6 6786410) |
| 158 | 2014 | 35238 | Appartement | 115000 | 79 | 1456 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350694.4 6786355) |
| 212 | 2014 | 35238 | Appartement | 109000 | 73 | 1493 | Rennes | 35238 | 352381204 | Les Cloteaux | POINT (350662.4 6786511) |
Agrégation
Il ne reste plus qu’à refaire l’agrégation en suivant la même procédure que pour les communes mais en retirant les iris où il y a eu moins de 5 ventes.
map_iris_dvf <- inter %>%
st_drop_geometry() %>%
group_by(CODE_IRIS, NOM_IRIS) %>%
filter(type == "Maison") %>%
summarise(nb = n(),
med_prixm2 = median(prixm2)) %>%
filter(nb >=5) %>%
right_join(map_iris_ctr) %>%
st_as_sf()
kable(head(map_iris_dvf,3))| CODE_IRIS | NOM_IRIS | nb | med_prixm2 | NOM_COM | INSEE_COM | geometry |
|---|---|---|---|---|---|---|
| 352380101 | Cathédrale | 5 | 2680 | Rennes | 35238 | MULTIPOLYGON (((351427 6789… |
| 352380102 | Hoche | 15 | 4671 | Rennes | 35238 | MULTIPOLYGON (((352032.6 67… |
| 352380109 | Vieux Saint-Étienne | 17 | 5263 | Rennes | 35238 | MULTIPOLYGON (((351501.3 67… |
Cartographie
Il ne reste plus qu’à faire la carte en reprenant le programme utilisé pour les communes :
mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_iris_dvf, type = "choro",
var = "med_prixm2",
breaks=mybreaks, pal = mypal,
border="gray", lwd=0.3,
leg_title = "prix en €/m2", leg_val_rnd = 0,
leg_pos = "left")
mf_map(map_iris_dvf, type = "prop",
var = "nb", inches = 0.08,
col="black",border="white",
leg_title = "nb. de ventes",
leg_pos = "topleft")
mf_layout(title = "Les ventes de maison de 2013 à 2020",
frame = TRUE, credits = "Sources : IGN et DVF",
arrow = F )MOD1 : Analyse sociologique / Chi2
On se propose dans un premier temps de modéliser des caractéristiques non-spatiales des individus à partir des données du RP 2019. On va considérer ici uniquement des attributs catégoriels ce qui nous conduira à formuler des hypothèses aboutissant à la création d’un tableau de contingence et la réalisation d’un test du chi-2.
Définir le sujet
Soit le sujet : Qui sont les nouveaux propriétaires ? Quel est leur niveau de qualification ?
Définir le statut d’occupation
Propriétaire ? Locataire ? Occupant à titre gratuit ?
Définir la notion de “qualification” ?
Le diplôme le plus élevé ? le nombre d’années d’étude ?
Définir la date
Année 2019 uniquement ? Résultats du RP 2019 (2017-2021) ?
Définir la population cible
Personnes installées récemment ? depuis quand ?
Définir la zone d’étude
Commune centre ? Agglomération ? Aire urbaine ? …
Formuler des hypothèses
Justes ou fausses, les hypothèses permettent de cadrer l’analyse. Elles peuvent faire intervenir plusieurs variables
Diplôme et propriété
Les nouveaux propriétaires sont plus souvent diplômés.
Âge et propriété
Les propriétaires sont plus âgés que les locataires.
Propriété et territoire
Les propriétaires sont plus nombreux en zone rurale
Propriété, âge, diplôme, territoire
Les jeunes ménages sont locataires en ville avant de devenir propriétaires dans les zones périurbaines ou rurales. Ils sont plus ou moins proches du centre selon leur niveau de qualification qui détermine à long terme leur revenu …
Charger les données
On charge ici les données ménages de l’aire urbaine.
FALSE COMMUNE ARM IRIS ACHL AEMM
FALSE 1: 22036 ZZZZZ ZZZZZZZZZ A12 1973
FALSE 2: 22036 ZZZZZ ZZZZZZZZZ B11 1992
Préparer l’analyse
Soit la relation entre statut d’occupation (Y) et diplôme le plus élevé du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés
Le choix de la population de référence est important. Ici on va sélectionner les ménages installés depuis moins de 3 ans
la sélection ou le regroupement des diplômes est également important car cela va influer sur les résultats du test.
la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage dans les zones urbaines
Sélection des individus et des variables
tab_sel<- tab_ind %>%
filter(as.numeric(ANEM) < 3) %>%
select(INSEE_COM, ANEM, DIPLM,STOCD, IPONDL)
knitr::kable(head(tab_sel,4))| INSEE_COM | ANEM | DIPLM | STOCD | IPONDL |
|---|---|---|---|---|
| 22036 | 2 | 15 | 10 | 1.034483 |
| 22036 | 2 | 13 | 23 | 1.034483 |
| 22036 | 2 | 13 | 21 | 1.034483 |
| 22036 | 1 | 13 | 10 | 1.034483 |
Recodage des modalités
On cherche le code des modalités dans le fichier des métadonnées
meta<-readRDS(paste0(myrep,"men_meta.RDS"))
metasel <- meta %>%
filter(COD_VAR %in% c("STOCD", "DIPLM"))
kable(metasel[,c(1,3,4)])| COD_VAR | COD_MOD | LIB_MOD |
|---|---|---|
| DIPLM | 01 | Pas de scolarité ou arrêt avant la fin du primaire |
| DIPLM | 02 | Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège |
| DIPLM | 03 | Aucun diplôme et scolarité jusqu’à la fin du collège ou au-delà |
| DIPLM | 11 | CEP (certificat d’études primaires) |
| DIPLM | 12 | BEPC, brevet élémentaire, brevet des collèges, DNB |
| DIPLM | 13 | CAP, BEP ou diplôme de niveau équivalent |
| DIPLM | 14 | Baccalauréat général ou technologique, brevet supérieur, capacité en droit, DAEU, ESEU |
| DIPLM | 15 | Baccalauréat professionnel, brevet professionnel, de technicien ou d’enseignement, diplôme équivalent |
| DIPLM | 16 | BTS, DUT, Deug, Deust, diplôme de la santé ou du social de niveau bac+2, diplôme équivalent |
| DIPLM | 17 | Licence, licence pro, maîtrise, diplôme équivalent de niveau bac+3 ou bac+4 |
| DIPLM | 18 | Master, DEA, DESS, diplôme grande école niveau bac+5, doctorat de santé |
| DIPLM | 19 | Doctorat de recherche (hors santé) |
| DIPLM | YY | Hors résidence principale |
| STOCD | 00 | Logement ordinaire inoccupé |
| STOCD | 10 | Propriétaire |
| STOCD | 21 | Locataire ou sous-locataire d’un logement loué vide non HLM |
| STOCD | 22 | Locataire ou sous-locataire d’un logement loué vide HLM |
| STOCD | 23 | Locataire ou sous-locataire d’un logement loué meublé ou d’une chambre d’hôtel |
| STOCD | 30 | Logé gratuitement |
On recode les modalités des deux variables en regroupant certaines.
tab_sel$STOCD<-as.factor(tab_sel$STOCD)
levels(tab_sel$STOCD)<-c("Proprétaire","Locataire",
"Locataire","Locataire",NA)
tab_sel$DIPLM<-as.factor(tab_sel$DIPLM)
levels(tab_sel$DIPLM) <- c("Aucun","Aucun","Aucun",
"BEP","BEP","BEP", "BAC","BAC",
"BAC+123","BAC+123","> BAC+3","> BAC+3",NA)
knitr::kable(head(tab_sel,3))| INSEE_COM | ANEM | DIPLM | STOCD | IPONDL |
|---|---|---|---|---|
| 22036 | 2 | BAC | Proprétaire | 1.034483 |
| 22036 | 2 | BEP | Locataire | 1.034483 |
| 22036 | 2 | BEP | Locataire | 1.034483 |
Tableau de contingence (FAUX)
Le plus simple semble être l’instruction table()
| Aucun | BEP | BAC | BAC+123 | > BAC+3 | Sum | |
|---|---|---|---|---|---|---|
| Proprétaire | 637 | 4132 | 3500 | 5738 | 3305 | 17312 |
| Locataire | 2944 | 9422 | 13158 | 12589 | 6237 | 44350 |
| Sum | 3581 | 13554 | 16658 | 18327 | 9542 | 61662 |
Tableau de contingence (JUSTE)
On pondère avec wtd.table() du package questionr.
programme
library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$STOCD,tab_sel$DIPLM,
weights = tab_sel$IPONDL)
knitr::kable(round(addmargins(tab_cont_wtd),0))| Aucun | BEP | BAC | BAC+123 | > BAC+3 | Sum | |
|---|---|---|---|---|---|---|
| Proprétaire | 751 | 4779 | 4167 | 7479 | 5125 | 22301 |
| Locataire | 4133 | 12437 | 21622 | 20839 | 11420 | 70450 |
| Sum | 4884 | 17216 | 25789 | 28317 | 16545 | 92751 |
Comparaison des résultats
- Tableau non pondéré … légèrement faux !
| Aucun | BEP | BAC | BAC+123 | > BAC+3 | All | |
|---|---|---|---|---|---|---|
| Proprétaire | 17.8 | 30.5 | 21 | 31.3 | 34.6 | 28.1 |
| Locataire | 82.2 | 69.5 | 79 | 68.7 | 65.4 | 71.9 |
| Total | 100.0 | 100.0 | 100 | 100.0 | 100.0 | 100.0 |
- Tableau pondéré … juste !
| Aucun | BEP | BAC | BAC+123 | > BAC+3 | All | |
|---|---|---|---|---|---|---|
| Proprétaire | 15.4 | 27.8 | 16.2 | 26.4 | 31 | 24 |
| Locataire | 84.6 | 72.2 | 83.8 | 73.6 | 69 | 76 |
| Total | 100.0 | 100.0 | 100.0 | 100.0 | 100 | 100 |
Visualisation du tableau de contingence
On choisit l’orientation du tableau et on l’affiche avec plot()
Tant qu’à faire, on améliore la figure avec des paramètres supplémentaires :
plot(mytable,
main = "Propriété et diplôme dans l'aire urbaine",
sub = "Source : INSEE - RP 2019 - Ménages installés depuis - de 3 ans",
col=c("lightyellow","lightgreen"))Test du Chi-deux
Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :
Pearson's Chi-squared test
data: mytable
X-squared = 1731.7, df = 4, p-value < 2.2e-16
- Interprétation : La relation est très significative (p<0.001)
Prolongements
On peut introduire toute une série de variantes à partir de ce premier programme. Par exemple, si on retire la commune-centre, les résultats sont assez différents :
tab_sel2<- tab_sel %>%
filter(INSEE_COM != codectr)
mytable<-wtd.table(tab_sel2$DIPLM,tab_sel2$STOCD,weights = tab_sel2$IPONDL)
plot(mytable,
main = "Propriété et diplôme hors commune-centre",
sub = "Source : INSEE - RP 2019 - Ménages installés depuis - de 3 ans",
col=c("lightyellow","lightgreen"))
Pearson's Chi-squared test
data: mytable
X-squared = 733.11, df = 4, p-value < 2.2e-16
A vous de deviner pourquoi :-)
MOD2 : Analyse spatiale / Gradient
On se propose d’analyser la variation d’un phénomène quantitatif (le prix de vente des maisons en € par m2) en fonction de la distance à un point remarquable (le centre de l’agglomération). On fait l’hypothèse que ce prix correspond à un gradient de décroissance, conformément aux théories de la rente foncière. Pour modéliser ce gradient, nous ferons appels à un modèle de régression linéaire ou non-linéaire.
Formuler l’hypothèse
Soit l’ensemble des ventes effectuées dans un rayon de 50 km autour d’une agglomération. Existe-t-il une relation entre (Y) le prix de vente mesuré en € par m2 de surface habitable et (X) la distance au point central de l’agglomération ?
N.B. Nous retirons la commune centre de l’analyse car on suppose que son comportement interne obéit à une logique différente
Définir l’espace d’étude
Nous allons retenir la distance à vol d’oiseau entre chaque point ou a eu lieu une vente et le centre de l’agglomération que l’on fixera en fonction de notre connaissance du terrain.
Calculer les distances au centre
On doit ajouter à notre fichier dvf la distance au point central retenu.
## Fichier des dvf
dvf<-readRDS(paste0(myrep,mydvf))
# Aire urbaine
map_com<-readRDS(paste0(myrep,"AA_mapcom.RDS"))
# Selection des dvf
dvf<-dvf %>% filter(INSEE_COM %in% map_com$INSEE_COM) %>%
filter(INSEE_COM != codectr)
# Transformation cartographique des dvf
map_dvf <- st_as_sf(dvf,coords =c("longitude","latitude"))
st_crs(map_dvf)<-4326
map_dvf<-st_transform(map_dvf,2154)
## Commune centre => centroïde
map_ctr <-readRDS(paste0(myrep,"CTR_mapcom.RDS"))
map_ctr <-st_centroid(map_ctr)
## Calcul de la distance
map_dvf$dist<-as.numeric(st_distance(map_dvf,map_ctr))/1000
## Visualisation rapide
plot(map_dvf[,"dist"])Modèle linéaire
On teste l’existence d’un modèle linéaire en fixant les paramètres suivants :
Paramètres
Ajustement
Call:
lm(formula = Y ~ X)
Residuals:
Min 1Q Median 3Q Max
-2309.1 -368.5 -47.0 314.5 5856.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2837.7774 5.9202 479.3 <2e-16 ***
X -45.1457 0.3075 -146.8 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 633.1 on 50618 degrees of freedom
Multiple R-squared: 0.2986, Adjusted R-squared: 0.2986
F-statistic: 2.155e+04 on 1 and 50618 DF, p-value: < 2.2e-16
Modèle exponentiel
On teste l’existence d’un modèle exponetielle en fixant les paramètres suivants :
Paramètres
Ajustement
Call:
lm(formula = Y ~ X)
Residuals:
Min 1Q Median 3Q Max
-1.95482 -0.14476 0.02937 0.19630 1.68254
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.9791861 0.0031466 2535.8 <2e-16 ***
X -0.0245287 0.0001635 -150.1 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3365 on 50618 degrees of freedom
Multiple R-squared: 0.3079, Adjusted R-squared: 0.3079
F-statistic: 2.252e+04 on 1 and 50618 DF, p-value: < 2.2e-16
Modèle logarithmique
On teste l’existence d’un modèle logarithmique en fixant les paramètres suivants :
Paramètres
Ajustement
Call:
lm(formula = Y ~ X)
Residuals:
Min 1Q Median 3Q Max
-2820.5 -364.0 -31.8 320.7 5900.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3801.054 12.449 305.3 <2e-16 ***
X -649.871 4.559 -142.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 638.5 on 50618 degrees of freedom
Multiple R-squared: 0.2864, Adjusted R-squared: 0.2864
F-statistic: 2.032e+04 on 1 and 50618 DF, p-value: < 2.2e-16
Modèle puissance
On teste l’existence d’un modèle puissance en fixant les paramètres suivants :
Paramètres
Ajustement
Call:
lm(formula = Y ~ X)
Residuals:
Min 1Q Median 3Q Max
-2.13041 -0.14325 0.04097 0.20474 1.55332
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.467618 0.006720 1260.1 <2e-16 ***
X -0.339949 0.002461 -138.1 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3447 on 50618 degrees of freedom
Multiple R-squared: 0.2738, Adjusted R-squared: 0.2738
F-statistic: 1.908e+04 on 1 and 50618 DF, p-value: < 2.2e-16
Choix du meilleur modèle
On décide ici de retenir le modèle 2 (exponentiel) et de stocker les valeurs estimées et résiduelles des prix au m2.
Agrégation par commune
On agrège les valeurs estimées et résiduelles par commune
mod_com <- st_drop_geometry(map_dvf) %>%
group_by(INSEE_COM) %>%
summarise(nb=n(),
prixm2 = mean(prixm2),
prixm2_estim = mean(prixm2_estim),
prixm2_resid = mean(prixm2_resid)) %>%
arrange(prixm2_resid)On en déduit les communes où le prix est plus bas que ce que laisserait prévoir la distance au centre :
| INSEE_COM | nb | prixm2 | prixm2_estim | prixm2_resid |
|---|---|---|---|---|
| 35281 | 1878 | 2435 | 2999 | -564 |
| 35261 | 18 | 1085 | 1571 | -486 |
| 35242 | 37 | 1045 | 1522 | -476 |
| 35244 | 31 | 1093 | 1535 | -443 |
| 35290 | 75 | 1122 | 1554 | -432 |
| 35164 | 80 | 1164 | 1577 | -413 |
Et celles où le prix est plus élevé que ce que laisserait prévoir la distance au centre
| INSEE_COM | nb | prixm2 | prixm2_estim | prixm2_resid |
|---|---|---|---|---|
| 35001 | 655 | 2401 | 2199 | 201 |
| 35069 | 1059 | 2225 | 2016 | 209 |
| 35177 | 486 | 2383 | 2162 | 222 |
| 35152 | 766 | 2200 | 1956 | 243 |
| 35047 | 2012 | 2533 | 2270 | 263 |
| 35051 | 1560 | 3134 | 2720 | 414 |
Cartographie des résultats
En se limitant à l’aire urbaine, on va visualiser les résultats du modèle. On commence par effectuer la jointure entre les résultats du modèle de régression et le fonds de carte des communes :
On peut alors cartographier les prix observés et théoriques :
mybreaks <-c(1000,1500, 1750, 2000, 2250, 2500, 2750, 3000, 6000)
mypal <- brewer.pal(n = 8, name="YlOrRd")
par(mfrow = c(1,2))
mf_map(map_com,
type = "choro",
var="prixm2_estim",
breaks=mybreaks,
pal=mypal,
leg_val_rnd = 0,
col_na = "gray80",
leg_pos = "topleft")
mf_layout("Prix théoriques",
arrow=F, frame=T,
credits = "Source : DVF, 2013-2020")
mf_map(map_com,
type = "choro",
var="prixm2",
breaks=mybreaks,
pal=mypal,
leg_val_rnd = 0,
col_na = "gray80",
leg_pos = "topleft")
mf_layout("Prix observés",
arrow=F, frame=T,
credits = "Source : DVF, 2013-2020")Et en déduire la carte des résidus :
mybreaks <-c(-700, -500, -300,-100,100,300,500,700)
mypal <- brewer.pal(n = 7, name="RdBu")
par(mfrow = c(1,1))
mf_map(map_com,
type = "choro",
var="prixm2_resid",
breaks=mybreaks,
pal=mypal,
col_na = "gray80",
leg_val_rnd = 0,
leg_pos = "topleft")
mf_layout("Prix résiduels",
arrow=F, frame=T,
credits = "Source : DVF, 2013-2020")